********************************************************************************
** Descriptive analysis "Mums Go Online: Is the Internet Changing the Demand 
** for Healthcare?" by Amaral-Garcia, Nardotto, Propper and Valletti.
********************************************************************************

** Table 1: Summary statistics 

global demoglistLSOA allpeople share20_44 full_time part_time white highskill4
global demoglistMicro twin anemia breech cardiaclung cordprol diabetes cervix hyper previousCS

**************************************************
* Program for 2001
cap program drop TeXcheckLSOA2001
program TeXcheckLSOA2001, rclass
	count if treated==1
	local Ntreat=r(N)
	count if treated==0
	local Ncont=r(N)
	file open myfile using "$ppath\Tables\BalanceCheckFirst.tex", write replace
	file write myfile "\begin{tabular}{lccccccc}" _n
	file write myfile "\hline\hline" _n
	file write myfile "\emph{A. LSOA Demographic} & \multicolumn{4}{c}{\textbf{Full sample}}& \multicolumn{3}{c}{\textbf{Matched sample}}\\" _n
	file write myfile "\emph{characteristics} & \multicolumn{4}{c}{\textbf{N=32,482}} & \multicolumn{3}{c}{\textbf{N=2,418}}\\" _n
	file write myfile "\cmidrule(l{4pt}r{4pt}){2-5}\cmidrule(l{4pt}r{4pt}){6-8}" _n
	file write myfile "\textbf{Variable}&Mean& Std.&Min&Max&\emph{Close}&\emph{Far}&P-value\\" _n

	file write myfile "\hline" _n
	file write myfile "\emph{Census 2001}\\" _n
	qui {
		foreach variab of var $demoglistLSOA {
			local labelvar : var label `variab'
			file write myfile %9s "`labelvar'"
			file write myfile " & "
			sum `variab'
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			file write myfile %5.2f (r(sd))
			file write myfile " & "
			file write myfile %5.2f (r(min))
			file write myfile " & "
			file write myfile %5.2f (r(max))
			file write myfile " & "
			sum `variab' if treated==1
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			sum `variab' if treated==0
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			ttest `variab', by(treated)
			scalar pval=r(p)
			file write myfile %5.3f (pval)
			file write myfile " \\"_n
		}
	}
	file write myfile "\hline" _n
	file close myfile
end

**************************************************
* Program for 2011
cap program drop TeXcheckLSOA2011
program TeXcheckLSOA2011, rclass
	count if treated==1
	local Ntreat=r(N)
	count if treated==0
	local Ncont=r(N)
	file open myfile using "$ppath\Tables\BalanceCheckSecond.tex", write replace

	file write myfile "\emph{Census 2011}\\" _n
	qui {
		foreach variab of var $demoglistLSOA {
			local labelvar : var label `variab'
			file write myfile %9s "`labelvar'"
			file write myfile " & "
			sum `variab'
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			file write myfile %5.2f (r(sd))
			file write myfile " & "
			file write myfile %5.2f (r(min))
			file write myfile " & "
			file write myfile %5.2f (r(max))
			file write myfile " & "
			sum `variab' if treated==1
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			sum `variab' if treated==0
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			ttest `variab', by(treated)
			scalar pval=r(p)
			file write myfile %5.3f (pval)
			file write myfile " \\"_n
		}
	}
	file write myfile " \\"_n
	file write myfile %9s "Distance LSOA - LE (km)"
	file write myfile " & "
	sum distance
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	file write myfile %5.2f (r(sd))
	file write myfile " & "
	file write myfile %5.2f (r(min))
	file write myfile " & "
	file write myfile %5.2f (r(max))
	file write myfile " & "
	sum distance if treated==1
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	sum distance if treated==0
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	ttest distance, by(treated)
	scalar pval=r(p)
	file write myfile %5.3f (pval)
	file close myfile
end

**************************************************
* Program the micro data on deliveries
cap program drop TeXcheckDeliveries
program TeXcheckDeliveries, rclass
	file open myfile using "$ppath\Tables\BalanceCheckThird.tex", write replace
	file write myfile "\\\hline" _n
	file write myfile "\multicolumn{4}{l}{\emph{B. Mothers' characteristics}}\\" _n
	file write myfile "\emph{and risk factors}&\multicolumn{4}{c}{\textbf{Full sample}}& \multicolumn{3}{c}{\textbf{Matched sample}}\\" _n
	file write myfile "& \multicolumn{4}{c}{\textbf{N=7,033,942}} & \multicolumn{3}{c}{\textbf{N=522,751}}\\" _n
	file write myfile "\cmidrule(l{4pt}r{4pt}){2-5}\cmidrule(l{4pt}r{4pt}){6-8}" _n
	file write myfile "\textbf{Variable}&Mean& Std.&Min&Max&\emph{Close}&\emph{Far}&P-value\\" _n
	file write myfile "\hline" _n
	* First I do age
	file write myfile %9s "Age"
	file write myfile " & "
	sum patage if patage>12 & patage<60
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	file write myfile %5.2f (r(sd))
	file write myfile " & "
	file write myfile %5.0f (r(min))
	file write myfile " & "
	file write myfile %5.0f (r(max))
	file write myfile " & "
	preserve
	use FinalData\delivery_main.dta, clear
	sum patage if treated==1
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	sum patage if treated==0 
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	ttest patage, by(treated)
	scalar pval=r(p)
	file write myfile %5.3f (pval)
	file write myfile " \\" _n
	restore
	* I continue with the remaining variables
	qui {
		foreach variab of var $demoglistMicro {
			* First, the full sample
			local labelvar : var label `variab'
			file write myfile %9s "`labelvar'"
			file write myfile " & "
			sum `variab'
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			file write myfile %5.2f (r(sd))
			file write myfile " & "
			file write myfile %5.0f (r(min))
			file write myfile " & "
			file write myfile %5.0f (r(max))
			file write myfile " & "
			* Second, the paired LSOAs
			preserve
			use FinalData\delivery_main.dta, clear
			sum `variab' if treated==1
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			sum `variab' if treated==0
			file write myfile %5.2f (r(mean))
			file write myfile " & "
			ttest `variab', by(treated)
			scalar pval=r(p)
			file write myfile %5.3f (pval)
			file write myfile " \\"_n
			restore
		}
	}
	*file write myfile "\hline" _n
	* Now Csections
	file write myfile " \\"_n
	file write myfile %9s "C-section in 2000 (\%)"
	file write myfile " & "
	sum csec if hesyear==2000
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	file write myfile %5.2f (r(sd))
	file write myfile " & "
	file write myfile %5.0f (r(min))
	file write myfile " & "
	file write myfile %5.0f (r(max))
	file write myfile " & "
	preserve
	use FinalData\delivery_main.dta, clear
	sum csec if treated==1 & hesyear==2000
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	sum csec if treated==0 & hesyear==2000
	file write myfile %5.2f (r(mean))
	file write myfile " & "
	ttest csec if hesyear==2000, by(treated)
	scalar pval=r(p)
	file write myfile %5.3f (pval)
	file write myfile " \\" _n
	file write myfile "\hline\hline" _n
	file write myfile "\end{tabular}" _n
	restore
	file close myfile
end

*********** END of programs for Table 1 ****************************


****************************************************
* First part of the Table: Census 2001 

use TempData\MatchedLSOAsStrict.dta, clear
keep lsoa_code_1 distance_1 lsoa_code_2 distance_2
gen idx=_n
reshape long lsoa_code_ distance_ , i(idx) j(treat)
gen treated=treat==1
drop treat
rename lsoa_code_ lsoa_code
rename distance_ distance

merge 1:1 lsoa_code using Demog\LSOADemogENGWALES_INCOME.dta  // All 2418 matched
gen nation=substr(msoa_code,1,1)
table nation
drop if nation=="W"
drop _merge

la var allpeople "Population"
gen share20_44=people20_29+people30_44
la var share20_44 "Age 20-44 (\%)"

gen white=whitebritishp+whiteirishp+whiteotherwhitep
replace white=100 if white>100
la var white "White (\%)"
gen highskill=largeemployersandhighermanageria+higherprofessionaloccupations
la var highskill "High skill workers (\%)"
gen highskill2=professionaloccupations+skilledtradesoccupations+associateprofessionalandtechnica
la var highskill2 "High skill workers (\%)"
gen highskill3=largeemployersandhighermanageria+higherprofessionaloccupations+intermediateoccupations
la var highskill3 "High skill workers (\%)"
gen highskill4=financialintermediation+realestaterentingandbusinessacti+ ///
publicadministrationanddefence+education+healthandsocialwork
la var highskill4 "High skill workers (\%)"

rename meanageofpopulationinthearea meanage

la var part_time "Part-time workers (\%)"
la var full_time "Full-time workers (\%)"

rename FamilyMarriedEchildren fam_child
la var fam_child "Families with children"

* Run the program to produce the table for the paper
TeXcheckLSOA2001


****************************************************
* Second part of the Table: Census 2011

use TempData\MatchedLSOAsStrict.dta, clear
keep lsoa_code_1 distance_1 lsoa_code_2 distance_2
gen idx=_n
reshape long lsoa_code_ distance_ , i(idx) j(treat)
gen treated=treat==1
drop treat
rename lsoa_code_ lsoa_code
rename distance_ distance

merge 1:1 lsoa_code using Demog\LSOADemogEW
gen nation=substr(lsoa_code,1,1)
table nation
drop if nation=="W"
drop _merge

la var allpeople "Population"

gen share20_44=people20_29+people30_44
la var share20_44 "Age 20-44 (\%)"
replace white=100 if white>100

la var white "White (\%)"
gen highskill=largeemployersandhighermanageria+higherprofessionaloccupations
la var highskill "High skill workers (\%)"
gen highskill2=professionaloccupations+skilledtradesoccupations+associateprofessionalandtechnica
la var highskill2 "High skill workers (\%)"
gen highskill3=largeemployersandhighermanageria+higherprofessionaloccupations+intermediateoccupations
la var highskill3 "High skill workers (\%)"
gen highskill4=informationandcommunication+financialintermediation+ ///
realestaterentingandbusinessacti+professionalscientificandtechni+education+ ///
healthandsocialwork+administrativeandsupportservice
la var highskill4 "High skill workers (\%)"

la var part_time "Part-time workers (\%)"
la var full_time "Full-time workers (\%)"

TeXcheckLSOA2011

****************************************************
* Third part of the Table: Mothers' characteristics
forval Y=2000/2011 {
	use Discharges\DatasetDeliveries_`Y'.dta, clear
	keep if epiorder==1
	keep patage csec del_csecelect previousCS twin anemia breech cardiaclung ///
	del_csecemerg cordprol diabetes cervix hyper hesyear soal
	foreach vari in csec del_csecelect del_csecemerg previousCS twin anemia breech cardiaclung cordprol diabetes cervix hyper {
		replace `vari'=`vari'*100
	}
	saveold TempData/TempDesc`Y'.dta, replace
}
use TempData/TempDesc2000.dta, clear
forval Y=2001/2011 {
	append using TempData/TempDesc`Y'.dta
}

la var patage "Age"
la var twin "Twins (\%)"
la var anemia "Anemia (\%)"
la var breech "Breech position (\%)"
la var cardiaclung "Cardiac and lung (\%)"
la var cordprol "Cord (\%)"
la var cervix "Cervix (\%)"
la var diabetes "Diabetes (\%)"
la var hyper "Hypertens. eclampsia (\%)"
la var previousCS "Previous C-section (\%)"
la var csec "C-sec in 2000 (\%)"

* Run the program
TeXcheckDeliveries

* Extra-check: number of deliveries
preserve
bysort soal hesyear: gen totmoms=_N
bysort soal: keep if _n==1
sum totmoms
restore

use FinalData\delivery_main.dta, clear
bysort id hesyear: gen totmoms=_N
bysort id: keep if _n==1
bysort id hesyear: keep if _n==1
ranksum totmoms if hesyear==2000, by(treated)
ranksum totmoms if hesyear==2011, by(treated)

/* Notice that the data file with the deliveries of the paired LSOAs is loaded
within the program (in the preserve-restore part), so the baseline data is the
full sample. Make sure it is that data which is loaded when you run the program.
*/

* Table 1 ends here ++++++++++++

********************************************************************************
* Figure 1: C-sec rate over time

use TempData/TempDesc2000.dta, clear
forval Y=2001/2011 {
	append using TempData/TempDesc`Y'.dta
}
collapse csec del_csecelect del_csecemerg, by(hesyear)
la var csec "All C-sec"
la var del_csecelect "Elective C-sec"
la var emerg "Emergency C-sec"
la var hesyear "Year"

#delimit ;
graph twoway connected csec del_csecemerg del_csecelect hesyear, color(black g6 g10) 
lpattern(solid dash dot) title() msymbol(circle square triangle) 
ytitle("C-section rates (%)") legend(off) xscale(range(1999 2012))
graphregion(color(white)) ysize(2.5) xsize(4);
#delimit cr
graph export $ppath/Graphs/CsecOverYearNoLegend.png, width(1000) replace  // 4000
* Figure 1 ends here ++++++++++++


********************************************************************************
* Figure 2 and 3. 
* Fig 2. Panel A: Broadband internet penetration between 2000 and 2011
* Fig 2. Panel B: Average actual speed in Mbit/s between 2000 and 2011
* Fig 3. Panel A: QGIS elaboration
* Fig 3. Panel B: Average actual speed in Mbit/s between 2000 and 2011 close vs. far LSOAs


* Panel A
use $extradatapath/BBpenetration.dta, clear

#delimit ;
graph twoway connected bbpen year, color(black) 
lpattern(solid)
ytitle("Internet penetration (%)") xscale(range(1999 2012))
yscale(range(0 90)) graphregion(color(white)) ysize(2) xsize(3);
#delimit cr
graph export $ppath\Graphs\BBpenOverYearFullRange.png, width(1000) replace  // 4000

* Panel B
use TempData\MatchedLSOAsStrict.dta, clear

expand 12
sort newid

by newid: gen hesyear=_n
replace hesyear=1999+hesyear

* Now I merge (each group)
rename le_code_1 LE
merge m:1 LE using "extradatapath/ActivationPeriodsFull.dta"
keep if _merge==3
drop _merge

foreach variab of varlist llu_year-ADSL2plus_year {
	rename `variab' `variab'_1
}

rename LE le_code_1
rename le_code_2 LE
merge m:1 LE using "extradatapath/ActivationPeriodsFull.dta"
keep if _merge==3
drop _merge

foreach variab of varlist llu_year-ADSL2plus_year {
	rename `variab' `variab'_2
}

gen adsl_tier_1=0
replace adsl_tier_1=1 if hesyear>=ADSL_active_year_1
replace adsl_tier_1=2 if hesyear>=llu_year_1
replace adsl_tier_1=3 if hesyear>=ADSL2plus_year_1

gen adsl_tier_2=0
replace adsl_tier_2=1 if hesyear>=ADSL_active_year_2
replace adsl_tier_2=2 if hesyear>=llu_year_2
replace adsl_tier_2=3 if hesyear>=ADSL2plus_year_2

* I use the coefficients in Table 1, pag. 18 of Ahlfeldt, Koutroumpis and Valletti
gen speed_1=7.869-0.293*distance_1^2+0.058*distance_1^3-0.003*distance_1^4 if adsl_tier_1==1
replace speed_1=8.214-0.287*distance_1^2+0.07*distance_1^3-0.005*distance_1^4 if adsl_tier_1==2
replace speed_1=8.672-0.491*distance_1^2+0.141*distance_1^3-0.011*distance_1^4 if adsl_tier_1==3

gen speed_2=7.869-0.293*distance_2^2+0.058*distance_2^3-0.003*distance_2^4 if adsl_tier_2==1
replace speed_2=8.214-0.287*distance_2^2+0.07*distance_2^3-0.005*distance_2^4 if adsl_tier_2==2
replace speed_2=8.672-0.491*distance_2^2+0.141*distance_2^3-0.011*distance_2^4 if adsl_tier_2==3

replace speed_1=exp(speed_1)
replace speed_2=exp(speed_2)
replace speed_1=128 if adsl_tier_1==0 | speed_1<128
replace speed_2=128 if adsl_tier_2==0 | speed_2<128
replace speed_1=speed_1/1000
replace speed_2=speed_2/1000

rename hesyear year
merge m:1 year using $ppath\Graphs\BBpenetration.dta
replace bbpen=6 if year==2000
replace bbpen=6 if year==2001
replace bbpen=9 if year==2002

replace speed_1=speed_1*bbpen/100
replace speed_2=speed_2*bbpen/100

* I do the graphs

* Figure 2, panel B
preserve
collapse speed_1 speed_2, by(year)
egen avgspeed=rowmean(speed_1 speed_2 )
#delimit ;
twoway connected avgspeed year, color(black) lpattern(solid)
ytitle("Internet speed over time (Mbit/s)") xscale(range(1999 2012))
graphregion(color(white)) ysize(2) xsize(3);
#delimit cr
graph export $ppath/Graphs/InternetSpeedOverTimeAll.png, width(1000) replace  // 4000
restore

* Figure 3, panel B (split between the two grops)
preserve
collapse speed_1 speed_2, by(year)
egen avgspeed=rowmean(speed_1 speed_2 )
#delimit ;
twoway connected avgspeed year, color(black) lpattern(solid)
ytitle("Internet speed over time (Mbit/s)") xscale(range(1999 2012))
graphregion(color(white)) ysize(2) xsize(3);
#delimit cr
graph export $ppath/Graphs/InternetSpeedOverTimeAll.png, width(1000) replace  // 4000
restore


********************************************************************************
* Figure 4: C-section rates over time
use FinalData\delivery_main.dta, clear

preserve
collapse csec (sd) sdcsec=csec, by(treated periodh)
rename periodh hesyear
gen hesyear2=hesyear^2
gen hesyear3=hesyear^3 
gen hesyear4=hesyear^4
gen hesyear5=hesyear^5

reg csec hesyear hesyear2 hesyear3 hesyear4 hesyear5 if treated==1
predict csechatreat

reg csec hesyear hesyear2 hesyear3 hesyear4 hesyear5 if treated==0
predict csechatcont

gen csechat=csechatreat if treated==1
replace csechat=csechatcont if treated==0

sort hesyear treated
by hesyear: gen dcsec=csec-csec[_n-1]
twoway (lpoly dcsec hesyear , bw(2) lc(gs6) lw(thick) lp(shortdash) ) , legen(off) ///
ytitle("C-section rate (%)")  ///
xlabel(7 "2000" 11 "2002" 15 "2004" 19 "2006" ///
23 "2008" 27 "2010" 31 "2012") xtitle("Financial Year") ///
graphregion(color(white)) 

twoway (lpoly csec hesyear if treated==0, bw(2) lc(gs6) lw(thick) lp(shortdash)) ///
(lpoly csec hesyear if treated==1, bw(2) lc(gs3) lw(thick)) /// 
(scatter csec hesyear if treated==0, mcolor(gs6) msymbol(D)) ///
(scatter csec hesyear if treated==1, mcolor(gs3)), legen(off) ///
ytitle("C-section rate (%)") xscale(r(6 31)) yscale(r(20 25.4)) ///
xlabel(7 "2000" 11 "2002" 15 "2004" 19 "2006" ///
23 "2008" 27 "2010" 31 "2012") xtitle("Financial Year") ///
graphregion(color(white)) 
graph export $ppath/Graphs\CsecTreatedControls.png, width(1000) replace

restore

********************************************************************************
* Figure 5: C-section rates over time (Split by high-low Education/Income)

use FinalData\delivery_main.dta, clear

gen elec=del_csecelec
la var elec "Elective C-sec"

gen emerg=del_csecemerg
la var emerg "Emergency C-sec"

egen hospnr=group(procode) // We need it for hospitals' FEs

gen experienced=prevpreg 
la var experienced "Multiple-time"

gen expdid=experienced*did
la var expdid "Multiple $\times$ Did"

gen expd2006=experienced*d2006

#delimit ;
global cov1="
agedum2 agedum3 agedum4 agedum5 agedum6 anemia breech cardiaclung cordprol 
diabetes disprop distress herpes hypereclamp cervix obesity  
placprevia prematmem previousCS prevabortion renal rh twins white";
#delimit cr

preserve
collapse csec (count) ngroup=csec, by(hesyear high_educmed)
gen csec1=csec/100
gen hi_bound=csec+(1.96*sqrt(csec1*(1-csec1)/ngroup))*100
gen low_bound=csec-(1.96*sqrt(csec1*(1-csec1)/ngroup))*100

#delimit ;
twoway 
(connected csec hesyear if high_educmed==1, color(black) lpattern(solid) )
(connected csec hesyear if high_educmed==0, color(black) lpattern(dash) )
(rcap low_bound hi_bound hesyear if high_educmed==1, lc(black) yline(0, lp(dash) lc(black)))
(rcap low_bound hi_bound hesyear if high_educmed==0, lc(black) yline(0, lp(dash) lc(black)))
,
title("C-section rates over time")
subtitle("Split by high-low education") 
ytitle("C-section rates (%)") xscale(range(1999 2012))
graphregion(color(white)) ysize(2) xsize(3) 
legend(on order(1 2) lab(1 "High education") lab(2 "Low education") lab(3 "") lab(4 ""));
#delimit cr
graph export $ppath/Graphs\CsecOverYearHighLowEducationCI.png, width(1000) replace  // 4000
restore

preserve
collapse csec (count) ngroup=csec, by(hesyear high_incmed)
gen csec1=csec/100
gen hi_bound=csec+(1.96*sqrt(csec1*(1-csec1)/ngroup))*100
gen low_bound=csec-(1.96*sqrt(csec1*(1-csec1)/ngroup))*100

#delimit ;
twoway 
(connected csec hesyear if high_incmed==1, color(black) lpattern(solid) )
(connected csec hesyear if high_incmed==0, color(black) lpattern(dash) )
(rcap low_bound hi_bound hesyear if high_incmed==1, lc(black) yline(0, lp(dash) lc(black)))
(rcap low_bound hi_bound hesyear if high_incmed==0, lc(black) yline(0, lp(dash) lc(black)))
,
title("C-section rates over time")
subtitle("Split by high-low income") 
ytitle("C-section rates (%)") xscale(range(1999 2012))
graphregion(color(white)) ysize(2) xsize(3) 
legend(on order(1 2) lab(1 "High income") lab(2 "Low income"));
#delimit cr
graph export $ppath/Graphs/CsecOverYearHighLowIncomeCI.png, width(1000) replace  // 4000
restore

restore

* Figure 5 ends here! ++++++++++++++++++++++++++++++++++++++++++++++++++++++


********************************************************************************
* Figure 6: From Redshaw, M. and Heikkila, K. (2010)


********************************************************************************
* Figure 7: Event Study
* Version 1b - csec rate
use FinalData\delivery_main.dta, clear

* Set as cov1 the following list of vars:
#delimit ;
global cov1="
agedum2 agedum3 agedum4 agedum5 agedum6 anemia breech cardiaclung cordprol 
diabetes disprop distress herpes hypereclamp cervix obesity  
placprevia prematmem previousCS prevabortion renal rh twins white";
#delimit cr

* All mothers (left panel)
preserve
areg csec $cov1 imd04ed imd04i dist_hosp i.hospnr, a(newid)
predict csec_res, res

bysort treated periodh: egen avgcsec=mean(csec)
bysort treated periodh: egen avgres=mean(csec_res)
bysort treated: gen totn=_N

bysort periodh treated: keep if _n==1

sort periodh treated
by periodh: gen dcsec=avgcsec-avgcsec[_n-1]
by periodh: gen dcsec_res=avgres-avgres[_n-1]

cap drop x s se
lpoly dcsec periodh , bw(2) lc(black) lw(thick) generate(x s) se(se) level(90)
qui replace se=se/100  // I have to rescale it
gen hi_bound=s+(1.28*se)*100
gen low_bound=s-(1.28*se)*100

twoway (line s x ,  lc(black) lw(thick)) ///
(line low_bound x , lp(dash) lc(gs6) ) ///
(line hi_bound x , lp(dash) lc(gs6) ) , legen(off) ///
ytitle("Difference in C-section rate (%) - Close - Far LSOAs")  yscale(r(-0.4 0.1 0.6)) ///
xlabel(7 "2000" 11 "2002" 15 "2004" 19 "2006" ///
23 "2008" 27 "2010" 31 "2012") xtitle("Financial Year") ///
yline(0) graphregion(color(white)) 
graph export $ppath/Graphs/EventStudy.png , width(1000) replace  // 4000

restore

* First time mothers (right panel)
preserve
areg csec $cov1 imd04ed imd04i dist_hosp i.hospnr if experience==0 , a(newid)
predict csec_res, res

keep if experience==0

bysort treated periodh: egen avgcsec=mean(csec)
bysort treated periodh: egen avgres=mean(csec_res)
bysort treated: gen totn=_N


bysort periodh treated: keep if _n==1

sort periodh treated
by periodh: gen dcsec=avgcsec-avgcsec[_n-1]
by periodh: gen dcsec_res=avgres-avgres[_n-1]

cap drop x s se
lpoly dcsec periodh , bw(2) lc(black) lw(thick) generate(x s) se(se) level(90)
qui replace se=se/100  // I have to rescale it
gen hi_bound=s+(1.28*se)*100
gen low_bound=s-(1.28*se)*100

twoway (line s x ,  lc(black) lw(thick)) ///
(line low_bound x , lp(dash) lc(gs6) ) ///
(line hi_bound x , lp(dash) lc(gs6) ) , legen(off) ///
ytitle("Difference in C-section rate (%) - Close - Far LSOAs")  yscale(r(-0.4 0.1 0.6)) ///
xlabel(7 "2000" 11 "2002" 15 "2004" 19 "2006" ///
23 "2008" 27 "2010" 31 "2012") xtitle("Financial Year") ///
yline(0) graphregion(color(white)) 
graph export $ppath/Graphs/EventStudyFirstTime.png , width(1000) replace  // 4000

restore


***************************************************************
** Fraction of mothers same hospitals
use FinalData\delivery_main.dta, clear
preserve
bysort newid treated procode_id : keep if _n==1
keep newid treated procode_id
keep if treated==1
save $extradatapath/ListHospTreated.dta, replace
*saveold TempData/ListHospTreated.dta, replace
restore

preserve
bysort newid treated procode_id : keep if _n==1
keep newid treated procode_id
keep if treated==0
save $extradatapath/ListHospControls.dta, replace
*saveold TempData/ListHospControls.dta, replace
restore

* Compute the shares
use FinalData\delivery_main.dta, clear
count
local totmoms=r(N)
* Switch the treatment variable
replace treated=treated+1
replace treated=0 if treated==2
merge m:1 newid treated procode_id using $extradatapath/ListHospTreated.dta
drop if _merge==2
count if treated==1 & _merge==1
local nomatch1=r(N)
drop _merge
merge m:1 newid treated procode_id using $extradatapath/ListHospControls.dta
drop if _merge==2
count if treated==0 & _merge==1
local nomatch2=r(N)
drop _merge

local samehosp=1-(`nomatch1'+`nomatch2')/`totmoms'
disp "Fraction of mothers same hospital: `samehosp'"




